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Abstract 



We study the quantum phase transition between a band ("ionic") insu- 
lator and a Mott-Hubbard insulator, realized at a critical value U = U r in a 



o 

ON. 

0> , bipartite Hubbard model with two inequivalent sites, whose on-site energies 



differ by an offset A. The study is carried out both in D = 1 and D = 2 
(square and honeycomb lattices), using exact Lanczos diagonalization, finite- 



C . 

O , size scaling, and Berry's phase calculations of the polarization. The Born 

O 

effective charge jump from positive infinity to negative infinity previously dis- 
covered in D = 1 by Resta and Sorella is confirmed to be directly connected 
with the transition from the band insulator to the Mott insulating state, in 
agreement with recent work of Ortiz et al. In addition, symmetry is anal- 
ysed, and the transition is found to be associated with a reversal of inversion 
symmetry in the ground state, of magnetic origin. We also study the D = 1 
excitation spectrum by Lanczos diagonalization and finite-size scaling. Not 
only the spin gap closes at the transition, consistent with the magnetic nature 
of the Mott state, but also the charge gap closes, so that the intermediate state 
between the two insulators appears to be metallic. This finding, rationalized 
within unrestricted Hartree-Fock as due to a sign change of the effective on- 



site energy offset A for the minority spin electrons, underlines the profound 

difference between the two insulators. The band-to- Mott insulator transition 

is also studied and found in the same model in D = 2. There too we find 

an associated, although weaker, polarization anomaly, with some differences 

between square and honeycomb lattices. The honeycomb lattice, which does 

not possess an inversion symmetry is used to demonstrate the possibility of 

an inverted piezoelectric effect in this kind of ionic Mott insulator. 
75.10.Jm, 71.20.Ad, 71.27.+a 
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I. INTRODUCTION 



Standard discussions of the Mott-Hubbard transition are generally concerned with lat- 
tices of equivalent sites. At zero temperature, the metal-insulator transition develops when 
the on-site electron-electron repulsion U reaches some critical value U c which, usually, also 
corresponds to the onset of characteristic magnetic correlations. In this paper we are con- 
cerned with the less common case where the system is ionic, encompassing "anions" and 
"cations". Earlier workers, including Nagaosa and Takimoto |1, and Egami, Ishihara and 
Tachiki fl[] considered the simplest two-site ionic generalization of the Hubbard model (hence- 
forth dubbed Ionic Hubbard Model, IHM), which exhibits the transition from a band insu- 
lator to a Mott insulator. Such a system is, in the absence of electron-electron repulsions, U 
= 0, a regular band insulator . As U increases above some critical value U c , a band insulator 
- Mott-Hubbard insulator transition is expected to take place. When U is sufficiently large, 
the inequivalence between anion and cation should in fact become irrelevant, and the ground 
state of a large U system of equivalent sites is Mott-insulating, with antiferromagnetic cor- 
relations. However, owing to the residual inequivalence of the two ionic sites, it will also 
exhibit other properties, which have only partly been explored so far. Of special interest 
is the anomalous behaviour of the polarization of the ionic solid across the transition. A 
very interesting quantity in this regard is the Born effective charge Z* associated with an 
infinitesimal "dimerizing" displacement of the ionic lattice, corresponding to a frozen q = 
optical phonon. Resta and Sorella |§ studied the IHM in D = 1 using Lanczos diagonal- 
ization and found, by the Berry phase method, very striking indications of a polarization 
anomaly at U — U c . Ortiz, Ordejon, Martin and Chiappe proposed a simple Hartree-Fock 
explanation to the anomaly, namely that at a first order magnetic transition between an AF 
insulator and a band magnetic insulator, the polarization may also abruptly jump. 

A number of important questions are apparently still open at this stage, in particular 

1. what is the physical origin of the polarization anomaly at the band-Mott transition, 
and what is its connection with symmetry? 
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2. what is the nature precisely at U = U c , of the "threshold" state between the band and 
the Mott insulator? 

3. how does the polarization anomaly depend on the dimensionality? In particular, will 
it survive in 2D and 3D, and if so, with what strength? 

4. bearing in mind that most experiments measure just only Z* 2 , can we identify a simple 
experimentally accessible quantity which could signal the Z* anomaly in magnitude 
and sign? 

In this paper we set out to discuss principally these questions. We will do it by study- 
ing more closely the same simple ionic Hubbard model |l],|3],|4[] , in particular by analysing 
the symmetry of its possible ground states, by calculating effective charges and excitation 
spectra, and generally by seeking to understand its properties by exact diagonalizations 
supplemented when necessary by the simple Hartree-Fock approximation. 

Firstly, we find that the Mott insulator does, for appropriate boundary conditions, pos- 
sess odd symmetry under inversion, contrary to the ionic band insulator, which has even 
parity. Secondly, Hartree-Fock reveals that one clue to the polarization anomaly lies in an 
effective reversal of the on-site energy offset (accompanied with the vanishing of the associ- 
ated band gap), taking place for the minority spin species only at U c . Thirdly, in agreement 
with Hartree-Fock, exact results suggest that at U c not only the spin gap, but also the true 
charge gap vanishes , indicating a metallic "threshold" state poised precisely at the brink 
between the band insulator and the Mott insulator. Fourthly, a fresh study of two differ- 
ent 2D lattices, namely square and honeycomb, shows that a direct band insulator-Mott 
Hubbard insulator transition and the associated polarization anomaly may survive in higher 
dimensions too, particularly in the 2D honeycomb lattice. Here, the anomaly is weaker than 
in D=l ( a jump instead of a divergence). Finally, we propose the piezoelectric effect in a 
non-centrosymmetric lattice, here exemplified precisely by the 2D honeycomb lattice, as the 
experimentally accessible quantity that will directly and strikingly change sign at the band 
insulator-Mott insulator transition. 
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The polarization calculations are carried out using the Berry phase technique, first intro- 
duced by King-Smith and Vanderbilt |J and by Resta |J , and further extended to a general 
many-body case by Ortiz and Martin |7]], and applied to the IHM by Resta and Sorella. 

The rest of this paper is structured as follows. In section [TT] we introduce the ionic 
Hubbard model and describe briefly the classical limit, t = 0, useful for understanding later 
the full quantum case. In section [H| we first present a discussion of the magnetic behavior 
expected by the model, followed by a more detailed discussion of symmetry, and by a finite- 
size study of the level crossings connected with the transition at U c . Section [IV| contains the 
Hartree-Fock theory of the band- Mott insulator transition, and of the polarization anomaly. 
Section [V] presents the full many-body calculation of the polarization, done by means of 
Lanczos diagonalization plus Berry phase, here specialized to the 2D honeycomb case, and 
finally in Section [V| we give a discussion of the possible detectability of the transition via 
the piezoelectric effect, followed by our general conclusions. 

II. IONIC HUBBARD MODEL 

We consider the Hubbard hamiltonian in ID (linear chain) and 2D bipartite lattices 
(square and honeycomb lattices). All lattices being bipartite, they are composed of A and 
B sublattices. To simulate ionicity, A and B are made inequivalent by onsite energies y and 
— y respectively. Because of the energy difference A between the A and the B sublattices, 
the electron population of the A sublattice is less than that of the B sublattice for A > 
(or viceversa for A < 0). The sublattices are connected by electron hopping . We assume 
a filling of one electron per site with equal numbers of spin up and of spin down electrons. 
Sites of the A sublattice are denoted by Ra, those of the B sublattice by Rb = Ra + Cm; 
where the vectors £ M connect an A site with its neighbouring B sites, \i — 1, 2, . . . , v. In the 
linear chain, where the length of the unit cell is 2, /i = 0,1, £o,i = ±1- In the square lattice 
(square side of unit length), \i = 0,1,2,3 and £0,2 = ±(1,0), £1,3 = ±(0,1)- In the honeycomb 
lattice (unit length is the side of the hexagon), \x = 0,1,2 and £ = (—5, ^r), £1 = (1,0), £2 
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(1) 



where we have used standard notation, and U is the Hubbard onsite electron-electron re- 
pulsive interaction. Electrons hop with matrix elements — > 0) between neighbouring 
sites from the A to the B sublattice along the £ M directions. We denote the ground state 
energy of our iV-electron system as E(N). When in need to distinguish between the number 
of spin up and spin down electrons, the ground state energy will be denoted by E(N |, iV |). 
Unless otherwise specified iV f= iV |= N/2. We shall study the charge gap and the spin 
gap of the system defined respectively as: AE c h arg e = E(N + 1) + E(N — 1) — 2E(N) and 
AE spm = E(N T +1, N i -1) - E(N !, N I). 



A. Behavior in the classical limit, t = 

It is instructive at the outset to consider what happens when all = 0. As sketched 
in Fig(|l|) if we set the hopping = in the Hamiltonian (H) the model is classical. It 
has a first order transition at U — A with a change of the macroscopic polarization per 
unit cell AP = ea/2. Both the charge and spin gaps are finite and coincide in the region 
U < A where they are given by AE = A — U . For U > A the charge gap remains finite 
AE c harge = U — A while the spin gap vanishes. Note that precisely at U = A, where the 
transition takes place, both the charge and the spin gap vanish. This kind of transition 
s expected to persist for t > 0, where it takes a standard band insulator, with a charge 
and a spin gap, over to an antiferromagnetic insulator, with a charge gap, and gapless spin 
excitations for large U. 
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III. MAGNETISM, SYMMETRY, LEVEL CROSSING AND POLARIZATION 

JUMP 



In this section we discuss the polarization properties of an electron system close to an 
antiferromagnetic transition. As discussed by Ortiz and Martin antiferromagnetism appears 
to play a crucial role and remarkably affects the behaviour of the polarization. Q 

Let us consider a finite electron hopping t^. For U = 0, one electron per site, the model 
is in fact described by a completely filled, spin independent band, separated by a finite gap 
(~ A, for large A) from the second, empty band. In the opposite limit of large U the 
standard strong coupling analysis of the Hubbard model leads to a charge gap ~ U, and to 
the well known Heisenberg Hamiltonian || Hj for spins: 

Hj = J £ ■ Sj (2) 

<ij> 

with an antiferromagnetic superexchange coupling J = The mapping to the Heisenberg 
model implies that for large U the model has gapless spin excitations. In one dimension they 
have been named "spinons" ||, and can be derived from the exact Bethe ansatz solution 
of the ID Heisenberg model. || In higher dimensions, where antiferromagnetic long range 



order is believed to exist at T = in 2D [ 10—12 1 and 3D [13], the gapless modes are spin 



waves. 

A band insulator - Mott insulator transition should therefore occur at some finite coupling 
U c . Of course, U c will differ from its classical value U® = A. There will be in general 
the possibility of intermediate metallic phases covering a range of U values. Even if the 
insulator-insulator transition is direct, quantum fluctuations may drive its character, for 
example from first order to second order, or else may split it into more than one transition 
|18| . Quantum antiferromagnets with long range order for U > U c (or quasi long range order 
in one dimension) can be described at the critical level, by the well known non linear sigma 
model, with action 



S = ±fdx d+1 (Vn) 2 (3) 



where n is a unit vector describing the local orientation of the antiferromagnetic order 
parameter, and the coupling constant g depends on the microscopic parameters U, A of the 
model. This effective model is well known to have a second order transition predicting that 



a spin gap opens up continuously for g > g c , alias U <U C [p~2 ] . 

We will for simplicity assume in the following that the transition is unique and is always 
second order even though a Hartree-Fock calculation 0] in the one dimensional model and our 
own in the 2D square lattice indicates the opposite. In fact, the Hartree-Fock approximation 
may fail anyway to describe correctly the order of the transition, as it is not appropriate to 
describe the gapless spin-wave excitations of the model in the ordered phase. Contrary to 
the linear chain and the square lattice, the Hartree-Fock method however correctly predicts 



(see Sec. [TV]) a second order phase transition in the honeycomb lattice, at least for the 
parameter values studied here. 

The magnetic Mott insulator for U > U c has a charge gap, and is therefore fully described 
by the nonlinear sigma model, at least in more than one dimensions. This model does not 
present any further phase transition but the one at the critical coupling g c , which is thus 
related to U c . It is possible however that the charge transition to a band insulator might 
occur for U values different from U c , as suggested by Fabrizio et al |T8"|] . We shall return to 
this point below. 



A. D=l: symmetry, level crossing and critical U c 

In a finite system with -say- L sites a level crossing may occur for some particular 
boundary conditions, if allowed by symmetry. In particular, in the one dimensional model, 
it can be easily proved that there exists a finite value U C (L) at which the ground state 
undergoes a change in the eigenvalue of the inversion symmetry operator R. Inversion 
symmetry around the site i = is defined by the following relations: 

Rc la Rj = c L-l- i ,a for* = 0,l,..-L-l (4) 
R\0 > = |0 > (5) 
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where |0 > is the electron vacuum state, by definition invariant under inversion symmetry. 

The additional relation (|5]) is necessary to completely define the inversion transformation 
R in the whole Hilbert space. Inversion does not interchange the A and the B sublattice 
and clearly commutes with the hamiltonian ([]]) for any U, provided the boundary conditions 
are real , namely periodic or antiperiodic (see Sec. |Vp. 

Since R 2 = I, the identity, the inversion has obviously two eigenvalues, ±1. We will show 
in the following that the ground state \ipo >u satisfies R\ipo >u= \4>o >u for U = 0, whereas 
for large U there is a change of sign and R\ipo >u= — IV'Q >u, so that a level crossing must 
occur at an intermediate coupling U = U C (L). 

In the non-interacting system the ground state is a direct product of a spin up and a 
spin down Slater determinants, both possessing the same orbitals of the lowest band. Both 
Slater determinants have a definite parity R a = ±1 and the inversion eigenvalue of the 
global wavefunction is given by their product 

R = i? T x R l = 1 

. Hence the band insulating state is even under inversion. 

In the large-C/ Mott insulator instead we can use the mapping to the Heisenberg hamil- 
tonian (0) , whose ground state in arbitrary dimensions can be generally written as ||13|| : 



l^o >u^oo= Hh,i2,---iL/2)S h S h ■■■S iL/2 \F > 

where S~ = cjjC^ is the spin lowering operator at the site i, \F >= II c {i|0 > is the 
ferromagnetic state along the spin-up direction, and the wavefunction $ is real. $ is also 



subject to the well known ("Marshall sign") condition | 13fl , i.e. the sign of the wavefunction 



is determined by the number of spin flips in the B sublattice (i odd): 

L/2 

$(ii,i 2 ,---,^/2)(-l) fc=1 > 0. 

The action of the inversion symmetry R on the spin lowering operators can be easily found 
by applying the definition given in Eqs.(|],|D, namely RS^R^ = Sj^_ ^ _ ^ and thus R 



maps an element of the basis S^S^ ■ ■ ■ S ih/2 \F > to another one S if ■ ■ ■ S ti R\F > where 
i' = L — 1 — i and R\F >= Il c L-i-i|0 >= ±\F >. The overall sign ± in the latter 

i 

equation represents just the inversion eigenvalue of the ferromagnetic state \F > and can be 
determined using the canonical anticommutation rules to restore the order of the creation 
operators cj after the application of the inversion operator R to the ferromagnetic state \F >. 
Since inversion symmetry does not change the Marshall sign we arrive to the conclusion 
that the inversion eigenvalue of the Heisenberg wavefunction coincides with the inversion 
eigenvalue of the ferromagnetic state \F > which is simple to compute. 

Doing that we find that for U — ► oo, the inversion eigenvalue can change sign depending 
on the boundary conditions (b.c): 

R = (-1) L / 2 + 1 for periodicb.c. 

R = (—1)^/2 for antiperiodicb.c. (6) 

This finally proves our initial statement; in particular a level crossing from an even state 
to an odd one has to occur in a periodic ring with L = An or in an antiperiodic one with 
L = An + 2. On the other hand, there will not necessarily be a level crossing in a periodic 
ring with L = An + 2 and an antiperiodic one with L = An. In summary, we conclude that 
the demise of the band insulator occurs via a symmetry change, whose finite-size signature 
is a parity switch from even to odd in the appropriate boundary conditions. 

B. Consequences on the calculation of the polarization 

As will be discussed in Sec.([VD the change of polarization in a many-body system can 
be obtained using a form of averaging over the boundary conditions. Thus the averaging 
necessarily include both periodic (PBC) and antiperiodic (APBC) boundary conditions. 
Now, in one or in the other, depending on L, a level crossing will necessarily occur at some 
finite U C (L). ^From the theory of polarization a strong variation of the polarization can be 
expected as a function of U around U C (L) even in presence of a perturbations such as a 
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dimerization 5 (see Sec. (|V|)) - Therefore , within the hypothesis that there exists only a well 
defined Mott transition at a critical value U c in the thermodynamic limit, we may expect 
that: 

U C (L) -^U c for L -> oo (7) 

that is, the U value where the level crossing occurs for large size is just the critical value of 
the magnetic transition. 

C. Charge and Spin Gaps in D=l 

Understanding charge and spin excitation gaps is a crucial point. We have studied these 
quantities in the D=l function of U/t, performing calculations on finite rings with 

PBC or APBC. By considering the sequence of closed shells with one electron per site, there 
is no level crossing and a finite size scaling analysis can be safely applied to the charge and 
spin gaps (see the end of Sec.[TI] for their definitions). For a general finite size system, the 
lowest order correction to any gap should be of the form jp. We have used a three parameter 
fit 

Ai = ^ + ^ + | + ... (8) 

including also a higher order L~ 3 term, to improve the accuracy. In Fig. (0) we show 
the finite size calculations of the gaps for the 6 and the 12 site ring, as well as the result 
extrapolated to the thermodynamic limit with the finite size scaling (§). In Fig. (Q) we 
also present the results for the spin and charge gaps for open shell rings of 6 and 12 sites. 
The overall behavior of the gaps for the largest open shell ring (L=12) is in agreement with 
the infinite size extrapolation of the closed shell rings, supporting the validity of our finite 
size scaling. Starting with the band insulator at small U, and increasing U, we find that 
both charge and spin gaps are very close, and decrease together until they appear to close 
at some U c . For U > U c , the charge gap turns sharply upwards, while the spin gap doesn't. 
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Precisely at U c ~ 2.75t, our fit suggests finite size gap corrections of the form i which 
implies not only spin, but also charge gapless excitations. This being the case, the system 
at U c is metallic. The nature of this metal is unknown and deserves further investigation. 

/,From our calculations it is hard to say whether charge and spin gaps will vanish at 



exactly the same U c , or at two slighly different values, as very recently proposed in [|18 



Nonetheless it is suggestive that large finite size charge gaps and small spin gaps become 
slightly inverted after extrapolation, which goes precisely in the direction of a charge gap 
closing at a slightly smaller U than the spin gap. 

We defer all discussion of the possible two-transition scenario to the work of Fabrizio et 



al. Jl8| and we will not further dwell on it in this paper, where we consider for simplicity a 
single U c . 



D. Extension to higher dimensions: D=2 

In the previous analysis of a level crossing in the model ([I]) we did not explicitly use the 
exact Bethe ansatz solution of ID systems. In fact the result that the inversion symmetry R 
for large U has the same eigenvalue of the corresponding ferromagnetic state \F > remains 
valid also in D=2, an so does the evenness of the band insulator at small U. 

Unfortunately, unlike D=l, the D=2 inversion symmetry, transforming (x, y) —>■ (L — 
1 — x, L — 1 — y) leaves the ferromagnetic state invariant on a bipartite lattice. Thus, a 
level crossing cannot be argued based on identically the same symmetry argument. As it 
turns out, however, it is again possible to generalize the argument by using a more elaborate 
symmetry operator, which changes eigenvalue in going from the U = state to large U 
state. In the square lattice this symmetry operator is easily identified. It is the mirror 
symmetry across the diagonal of the square lattice L = I x / with 1/2 odd (see picture |]). 
In the honeycomb lattice one can also find a symmetry operation with the same property. 
However, it is much more involved and we will not discuss it in detail here. We only mention 
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that this symmetry is obtained by a 120° rotation around a site followed by an additional 
gauge transformation cj — > cje- 7 ^ with suitable angles t . [[14]] 

Based on this analysis, we can therefore conclude that, upon averaging over the boundary 
conditions, there will be, both in the square and in the honeycomb lattice, a level crossing 
when the system is in the Mott state, but none in the band insulator state . Therefore, we 
should expect a polarization anomaly, and a metallic state at the transition, also in these 
two dimensional cases. 

However, before moving on to do numerical work and check these expectations in these 
more difficult problems, it is wise to solve them in simple mean-field which, as the D=l case 
demonstrated, H is always very instructive. 



IV. 2D BIPARTITE LATTICES: HARTREE-FOCK APPROXIMATION 

We shall consider the Hubbard model on the bipartite honeycomb lattice, defined in Sec. 
(y), and on the simpler square lattice. In the latter case, in order to remove the nesting 
degeneracy of the non interacting 2D Fermi surface, we have also studied the effect of the 
next-nearest neighbor hopping. 

In the Hartree -Fock approximation the ground state of the Hamiltonian is approximated 
by a Slater determinant. We may further assume that this Slater determinant \SD > is 
factored in spin space 

\SD >= \SDi > ® \SD l > (9) 

With this choice it is simple to obtain the Hartree-Fock Hamiltonian by linearizing the 
interaction term: 

Ur^ni = U(< SD l \n l \SD l > n T + < ST> T |n T |ST> T > nfi-U < SD l \n l \SD l >< SD^\SD^ > 
In this way we obtain: 

H Hartree- Fock = + H{ + E const 
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where 

H a = e kc{ Ba c kA(T + h.c. + A a J2 ( c kAa c kAa - c{ Ba c kBa ) (10) 

keBZ k&BZ 

where we have employed a Fourier transform in the two sublattices and correspondingly we 
have defined the complex function: 

efc = £e^ (11) 
In the square lattice case the presence of the next-nearest hopping implies a further term: 

H t' = e 'k( c lAa c kAa + c\ B(T c kB(T ) (12) 

keBZ,cr 



with e' k = At' cos k x cos k y to be added to H a in Eq.(|lO|). For a uniform solution in both 
sublattices the average spin densities are given by: 

, pa i for Re A 
kSD^rISD^ ={ (13) 

p Bl for Re B 

. Pai for Re A 
< SD^ tR \SD^ > = { (14) 

p m for Re B 

The parameters Aj and A^, defining the Hartree -Fock hamiltonians are given therefore by: 

A T = y + j(Pai -Pbt) (15) 
A| = ^ + ^(PA 1 -PB 1 ) (16) 

The constant is obtained after little algebra 

E const = ^{pA\ + PB^)N 1 + Tf(pAl + pBl)N^ - ^~{pA\pA{ + pB\pBl) 

where iVj and N± are the total number of spin up and spin down particles (N^ = N± = L/2). 



At one electron per site all the bonding bands E kjCr = — J \e k \ 2 + A^ are occupied by the 
spin up and spin down electrons and the total Hartree-Fock energy is 

Etot = E const + £ E kja (17) 

keBZa 
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By minimizing the total energy we obtain two self consistent equations for the variational 
parameters A CT : 



A A_ CT „ 1 
A. = ¥ --F^E^ (18) 

1 L k&BZ ^k-a 



These equations can be easily solved by inserting a trial initial value for Aj and Aj in the 
rhs and iterating the result until selfconsistency is reached. For small U only a self consistent 
solution with Aj = Aj is possible, with a small charge transfer from the electron-rich site 
B to the electron-poor site A. The plot of the quantity pa — Pb = (Aj + Aj — A) 2 /U, is 
shown in Fig. (||). 

For large U a broken symmetry solution with Af ^ Aj is possible for U sufficiently large. 
A finite value of the staggered magnetization a = \{px\ ~ Pa\ ~ Pbi + Pbi) is stable, and is 
given by: 

A j. - A T = Ua 

For large enough U, A can be neglected and we approach asymptotically the standard 
antiferromagnetic solution, where 

At = -A ; . 

It is clear therefore that after the magnetic transition, which occurs first in HF, there is a 
slightly larger critical U when one of the two bands becomes gapless with A CT = 0. A plot of 
the Hartree-Fock band gaps Aj and Aj is shown in Fig. @. The self consistent equations 
(PI) for the 2D square lattice case are not modified by the presence of t' at least in the 
insulating HF phases where the bonding bands E ka = e' k — E^^ are full and the antibonding 
bands E$ a = e' k + E Ka are empty for both spin up and down electrons. However for t' 
we found a finite U interval - between the band insulator at small U and the Mott-Hubbard 
insulator at large U - where the bonding and the antibonding bands do overlap in a finite 
energy range, leading to a fully metallic behavior (see Fig.^). Here the above analysys should 
be slightly modified to take into account the gain in energy obtained by occupying some 
states of the antibonding band. 
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We will not study the square lattice tt' model further, since the general three-phase 
sequence: band insulator - extended metal phase - Mott-Hubbard insulator is not our main 
concern. And although the Hartree-Fock approximation is by no means exact, it certainly 
represents a good starting point in dimension larger than one, so that we expect the tt' 
model to have indeed an extended metal phase between the two insulators. 

According to Hartree-Fock, the two insulator phases are however still adjacent on the two 
sides of U c in the square lattice with first neighbor hopping (t'=0), and in the honeycomb 
lattice, with any hopping. In the former, the nesting property of the non interacting Fermi 
surface affects dramatically the HF solution, making it effectively one dimensional, and 
leading to an unphysical first order transition. |4| Moreover, even beyond Hartree-Fock, the 
first-neighbor square lattice model is not generic enough. An arbitrarily small next neighbor 
hopping will take it away to a tt', where the two insulators are most likely no longer adjacent. 

We conclude that, at least at the HF level, the simplest model which generically possesses 
a transition between a Mott and a band insulator in 2D is the honeycomb lattice. Here i) 
the HF transition to a magnetic phase is second order; ii) there is no nesting of the non 
interacting Fermi surface, which consists of just two k-points; and iii) the HF solution is 
always insulating (band or Mott) as a function of U with the exception of a point where 
there is semimetallic behavior. For the above reasons, rather well justified at the HF level, 
in 2D we shall restrict our study to the honeycomb lattice case. 

In the honeycomb HF solution, there is a small difference between the two critical values 
of U, one where magnetism sets in, and the next where a band gap closes, Fig. (||). Bearing 



in mind the discussion of Sec. [II [B|| , we cannot judge whether this difference and its sign 
is just an artifact of the HF calculation or a real one. Because of the crudeness of Hartree 
Fock, most likely the second. 

The evolution of the full Hartree-Fock band structure for increasing U (not shown) is also 
instructive. The bands display a full gap both in the band insulator, and in the magnetic 
insulator. The presence of a spin gap in the latter is clearly an artifact due to breaking 
of spin rotation invariance characteristic of Hartree-Fock. At the upper critical U, where 
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the charge gap closes, there is linear band dispersion, a point-like Fermi surface, and a 
Dirac massless spectrum centered at the 2D zone boundary. Therefore we could expect a 
semimetallic state to exist at the treshold between the two different insulators. 



A. Nature of the 2D polarization anomaly in the ID and 2D honeycomb lattice 

We argue that close to the point where the effective gap Aj changes sign, the polarization 
will change dramatically as a function of U. To verify our statement, we make use of 
the theory of the polarization, in the geometric phase formulation 0, ||, 0, |15| . In the 
Hartree-Fock case, this amounts to calculate the contribution to the polarization of the spin- 
polarized band which becomes gapless at the critical point, when one effective gap parameter 
changes sign. In order to better illustrate that, we consider a model system consisting of 
spinless noninteracting fermions, described by the Hamiltonian (|T0|)with a density of 1/2 
fermion per site. With periodic boundary conditions, we can calculate analytically the 
change of polarization as we vary continuously the fermion energy gap A, or alternatively 
as we introduce a small dimerizing distortion 5 in the hopping. In the linear chain, without 
any distortion, 5^0, the variation of the polarization as a function of the energy difference 
A between anions and cations is described by a step function : Pq(A) — |4 ( see Ref. f7j). 
When we vary A without changing its sign, the polarization does not change, while when 
the sign changes, the polarization jumps by 4f as anticipated. We introduce a dimerizing 
distortion ("optical phonon-like" ) 5 which amounts to a change of the hopping t — > t(l ± 5) 
for electrons hopping to the right and the left respectively of any B site. We can evaluate 
analytically the "Born effective charge", which is the derivative of the polarization with 
respect to 5 at 5 = : 

pci I 

limPi(A) = —{k'K ~^K + E)} (19) 

where k = (1 + ^iY 1/2 , k' = |(1 + ^)~ 1/2 . K = K{k) and E = E(k) are complete 
elliptic functions of the first and second kind. 
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In this model, we have thus recovered the divergent effective charge anomaly discovered 
numerically in the ID many-body calculation of Resta and Sorella ||. The same analysys 
holds for two dimensional or higher d-dimensional lattices provided that at the transition 
point A = the bonding and the antibonding bands touch at a d — 1 dimensional (Fermi-) 
surface, (nesting). That is the case in particular for the D = 2 square lattice but not for the 
honeycomb lattice. 

In the honeycomb lattice, using the geometric phase technique for the mentioned spinless 
fermion model, it is straightforward to verify that without any distortion the polarization 
does not change as we vary A. Introducing a distortion along the £i direction: to = ti = 
t(l — 8),tx = t(l + 25), one can then verify that the component of the change in polarization 
which is orthogonal to £i vanishes, at least in first order in S. 

Along £i, it is difficult to derive an analytic expression for the change in polarization . 
We can however extract the singular part of the effective charge, in the neighbourhood of 
the zero of A 

Bs P « A » = < 20 > 

A-»0 1 1 

The effective charge has therefore a finite symmetric jump when A changes sign. Since this 
is only the singular part, there will be in addition a smooth background contribution shifting 
this jump with respect to zero.( In D=l, the singularity was infinite, so there the shift was 
irrelevant). 

In conclusion, if the Hartree-Fock approximate picture were correct, we should expect 
such a polarization jump to emerge in the full many-body calculation for the 2D honeycomb 
lattice of the next section. As it will turn out, the jump is indeed verified, see Fig. 
Moreover, the quantitative size of the jump is remarkably close to twice the value we have 
just obtained analytically. This result may be understood by assuming that each of the 
two bands with opposite spins, are forced to close together their single particle gaps, due to 
the spin rotation invariance of the model. This also suggests that the jump in the effective 
charge in 2D, similarly to the jump in the polarization in ID, maybe due to topological 
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reasons and therefore rather general and weakly dependent on the parameters defining the 
model. 



V. 2D HONEYCOMB LATTICE: MANY-BODY CALCULATON OF THE 

POLARIZATION 

We now proceed to do proper polarization-or more precisely, effective charge- calculation 
for a 2D interacting case. We restrict our study of the Hamiltonian ([!]) to the honeycomb 
lattice ( a sort of 2D hexagonal BN), which is slightly more interesting because it does 
not possess inversion symmetry, and can have, e.g., a piezoelectric effect (see Sec.|V|). The 
sites of the A sublattice are given by Ra = R m ,n, with R m ,n = (| n + 3m, ^n)a, a is 



•■2 1 2 

the lattice constant and 1 < n < N y , 1 < m < Nx + 2 . The lattice is space periodic 



under translations 7\ = (|, ^)N y a and T 2 = (3, 0) Nx ^ 2 a. The reciprocal lattice vectors 
are G\ = ^=^(0,1) and G 2 = \/3). The sites of the B sublattice are given by: 

Rb = Rmn + A* = 0' 1) 2 where the three vectors ^ which connect neighouring A and 
B sites are £ = (-§, ^)a, Ci = (1, 0)a and f 2 = (-5, -^) a - Obviously Co + Ci + 6 = 0. 
We can also consider, for the purpose of mimicking a uniaxial stress along the £1 direction: 
to = t 2 = t(l — 5) and t% = t(l + 25). The Hamiltonian depends parametrically on 5, 
H = H(5). In order to calculate the difference in the polarization of the system as the 
parameter is varied between its values and 5, we consider families of Hamiltonians Hk(5) 
obtained from (HD by introducing a complex hopping and substituting — > t M e* fc '^, where 
the parameter k is given by k = k\G\ + k 2 G 2 , with < k a < 1, a = 1,2. This is 
equivalent to imposing generalised boundary conditions on the original Hamiltonian ||16|| : 



If ip(ri, . . . ,rj, . . .) is an eigenfunction of H, then generalized boundary conditions imply 
ip(ri, . . . , Tj + T a , • • •) = e t2nk °"ip(ri, . . . , Tj, . . .), a = 1,2. Periodic boundary conditions 
correspond to k a = , antiperiodic to k a = 1/2. The polarization difference between two 
states which are characterised by the initial and final value of the distortion, and S is given 
as an integral over the generalised boundary conditions: 



19 



nG a -AP = e [ dkp(T a (5)-T a (0)) 



Jo 



(21) 



where fl is the unit cell volume, a and (5 take alternatively the values 1 or 2, T a (S) is the 
many-body generalisation of the geometric phase: 



and $o(5, k) is the ground state of H k (S), which satisfies periodic boundary conditions. 

In the numerical calculation, we adopt a cell with 8 lattice sites (N x — N y — 2) and 8 
electrons. The number of spin-up and spin-down electrons is 4, and the dimension of the 
Hilbert space is only 4900. The system is really very small, but we expect that averaging 
over the boundary conditions will reduce the finite size effects. In the calculation we are 
interested to study the behaviour of AP, the polarization difference between the undistorted 
(5 = 0) and the distorted (5 ^ 0) case, for different values of the onsite interaction U . 

We use a discretized form to calculate the geometric phase numerically 



The geometric phase however is only defined modulo 2ir and an uncertainty in the result 
for the difference in polarization arises. There will be an ambiguity in the value of the 
polarization difference up to a quantum. However the polarization difference should remain 
unambiguous for neighboring values of 5. We can thus fix the 2n uncertainty in T(5) by 
requiring that F(8) is continuous in 5 for a given value of U. We may also require that it 
is continuous for neighboring values of U for fixed 5, so long as we do not cross a point of 
degeneracy in the electronic spectrum. 

Convergence in the integration with respect to kp has not been easy and we had to use grids 
with 150 or 300 points to get reliable results. At the end we found that in the undistorted 
case, similarly to the one-dimensional case, there is a critical point U c , where the geometric 
phase changes discontinuously by 7r for one kp in the 2D Brillouin Zone. As in one dimension 
this effect maybe related to the presence of a level crossing as a function of U/t for a particular 




(22) 




(23) 
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boundary condition. Introducing next the dimerizing distortion, we found that the effective 
charge, as in the one-dimensional case, is positive in the region below U c and negative above. 
At U c it is discontinuous, with a finite jump, instead of the one-dimensional divergence. A 
plot of the effective charge as a function of U/t is shown in Fig. (Q). 

VI. PIEZOELECTRIC INVERSION, AND CONCLUDING REMARKS 

The transition between an ionic band insulator and a Mott-Hubbard insulator as exhib- 
ited by the prototype model (Eq.0) is quite interesting. We have extended existing studies 
in D= 1, probing deeper into the transition. We have also carried out newer investigations 
in D=2, where same kind of transition is shown to take place, with some differences between 
the square and the honeycomb lattice. 

With reference to the list of questions presented in the introduction, we can now formulate 
the following answers. 

1. The physical origin of the polarization anomaly is connected with the symmetry switch 
between an even state, typical of the band insulator, to an odd state, typical of the Mott 
insulator. The symmetry switch, which we |§ and others [§]] had noted earlier for finite 
size, is established here for arbitrary size. Moreover, in D=l, even and odd refer to simple 
inversion symmetry, in D=2 the pertinent symmetry operation is different and depends on 
the lattice. 

2. The threshold state between the ionic and the Mott insulator, where the polarization 
abruptly changes sign, is one where the charge gap seems to close, and is therefore metallic. 
This result underlines the profound difference between the two types of insulator, and invites 



further studies, which are now beginning to appear |18J of this transition. 

3. In D=2 the same model has again a band-to-Mott insulator transition. The polariza- 
tion anomaly is generally weaker, to a degree which depends on the lattice. In one dimension 
the effective charge diverges at the transition undergoing an abrupt sign change from +oo to 
— oo. In two dimensions this sign change persists. We have shown that in the 2D honeycomb 



21 



lattice the divergence turns to a jump, again with a change of sign. Here the magnitude of 
this jump is found to coincide almost exactly with ea/n, suggesting that this jump could be 
an experimentally detectable quantity dependent only on the bulk lattice constant and no 
other details of the actual material. 

4. The total piezoelectric coefficient 7 might be used to detect the polarization anomaly, 
because it is sensitive to the sign of Z* . It is conventionally divided into two contribu- 
tions: the first, 70, is purely electronic, and is related to a uniform stress of the solid (a 
simple scaling of all distances in the unit cell); the second, more important contribution, 
is obtained by keeping the unit cell fixed and by changing the distance between atoms. In 
lattices without inversion symmetry, the piezoelectric coefficient 7 is directly proportional, 
in magnitude and sign, to the effective charge Z* , in the form (we omit here for simplicity 
tensorial and vectorial indices) |T9| , |2~0 



7 = 70 + Z*i (24) 
where the constant £ represents the internal strain parameter. 

An anomaly in Z* will clearly reflect in the piezoelectric coefficient. Piezoelectric mea- 
surements could therefore permit the detection of a band-to-Mott insulator transition, pro- 
vided the system does not possess inversion. The honeycomb lattice which we have studied 
in D=2 is the simplest one that does not possess inversion symmetry. Our study shows that 
the transition is not suppressed by piezoelectric strain. 

We have not been able to identify, at present, a likely compound where this kind of 
transition could be experimentally detected. It would seem possible, nevertheless, that 
suitable systems could be engineered, especially in the organic world. 
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FIGURES 

FIG. 1. The model at the classical level when the hopping t is set to zero. The energy 
difference between an A site and a B site is A. At 1 electron per site filling, when the onsite 
interaction U < A the B site is doubly occupied and the A site empty. When U > A both sites 
are single occupied. 

FIG. 2. Calculation on 12 site rings (squares) and 6 site rings (triangles) of the charge (black) 
and spin (white) gaps, using closed shells at fillings of 1 electron per site. The data of the gaps for 
the 8, 10 and 12 sites are used for a finite size scaling extrapolation to the infinite number of sites. 
The curves of the charge (solid line) and spin (dashed line) gaps are shown. The two curves seem 
to coincide in the region below U c 

FIG. 3. Calculation on 12 site rings (solid and long dashed line) and 6 site rings (short dashed 
and dotted lines) of the charge (long dashed and dotted lines) and spin (solid and short dashed 
lines) gaps, using open shells at fillings of 1 electron per site. The results of the 12 site calculation 
practically coincide with the finite size scaling results of the closed shells. 

FIG. 4. Square lattice of side I with 1/2 odd. The number of lattice sites which lie on either 
side of the diagonal is odd. 

FIG. 5. honeycomb lattice, Hartree-Fock calculation of the staggered magnetization a and the 
density difference pa — Pb- The transition is 2nd order. The arrow indicates the transition point. 

FIG. 6. honeycomb lattice, Hartree-Fock calculation of the effective energy gap A a . Two 
different critical values U&, U C 2 of the interaction are identified. Below U& there is only a solution 
with Aj = Aj_. Above U c \ a solution with Aj 7^ Aj exists leading to a finite value of the magne- 
tization a. One of the spin bands becomes gapless at a slightly different value of the interaction 

u c2 . 

FIG. 7. Square lattice with nearest neighbor hopping: HF calculation of the staggered magne- 
tization a, the density difference pa — Pb and the band insulator gap Eq 
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FIG. 8. honeycomb lattice, plot of the effective charge Z* as a function the Hubbard onsite 
interaction U for A = 1. The effective charge shows a finite jump at U c unlike the ID case where 
the jump diverged. 
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